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Mass- Conserved  Phase  Field  Models  for  Binary  Fluids 

Jie  Sheii!  Xiaofeng  YangTand  Qi  Wang* 


Abstract 

The  commonly  used  incompressible  phase  field  models  for  non-reactive,  binary  fluids,  in  which 
the  Cahn-Hilliard  equation  is  used  for  the  transport  of  phase  variables,  conserve  the  total  volume  of 
each  phase  as  well  as  the  material  volume,  but  do  not  conserve  the  mass  of  the  fluid  mixture  when 
the  densities  of  two  components  are  different.  In  this  paper,  we  formulate  the  phase  field  theory  for 
mixtures  of  two  incompressible  fluids,  consistent  with  the  quasi-compressible  theory  [28],  to  ensure 
the  conservation  of  mass  and  momentum  for  the  fluid  mixture  as  well  as  the  volume  for  each  fluid 
phase.  In  this  formulation,  the  mass-average  velocity  is  no  longer  divergence-free  (solenoidal)  when 
the  densities  of  two  components  in  the  mixture  are  not  equal,  making  it  a  compressible  model  subject 
to  an  internal  constraint.  An  efficient  numerical  method  is  then  devised  to  enforce  this  compressible 
internal  constraint.  Numerical  simulations  in  confined  geometries  for  both  the  compressible  and  the 
incompressible  models  are  carried  out  using  spatially  high  order  spectral  methods  to  contrast  the  model 
predictions.  Numerical  comparisons  show  that  (a)  the  predictions  by  the  two  models  agree  qualitatively 
in  the  situation  where  the  interfacial  mixing  layer  is  thin;  and  (b)  the  predictions  differ  significantly  in 
binary  fluid  mixtures  undergoing  mixing  with  a  large  mixing  zone.  The  numerical  study  delineates  the 
limitation  of  the  commonly  used  incompressible  phase  field  model  and  thereby  cautions  its  predictive 
value  in  simulating  well-mixed  binary  fluids. 


1  Introduction 

Phase  field  models  have  been  used  successfully  to  study  a  variety  of  interfacial  phenomena  like  equi¬ 
librium  shapes  of  vesicle  membranes  [12,  13,  14,  15,  16,  35],  blends  of  polymeric  liquids  [36,  37,  38,  17], 
multiphase  fluid  flows  [19,  23,  24,  28,  25,  41,  40,  42,  43,  44,  45],  dentritic  growth  in  solidification,  mi¬ 
crostructure  evolution  [21,  29,  22],  grain  growth  [8],  crack  propagation  [9],  morphological  pattern  forma¬ 
tion  in  thin  films  and  on  surfaces  [26,  30],  self-assembly  dynamics  of  two-phase  monolayer  on  an  elastic 
substrate  [27],  a  wide  variety  of  diffusive  and  diffusion-less  solid-state  phase  transitions  [10,  39],  dislo¬ 
cation  modeling  in  microstructure,  electro-migration  and  multiscale  modeling  [34].  Multiple  phase-field 
methods  can  be  devised  to  study  multiphase  materials  [40].  Recently,  phase  field  models  are  applied  to 
study  liquid  crystal  drop  deformation  in  another  fluid,  liquid  films,  polymer  nanocomposites,  and  biofilms 
[19,  23,  24,  28,  25,  41,  40,  42,  43,  44,  18,  46,  5], 

Comparing  to  other  mathematical  and  computational  technologies  available  for  studying  multi-phase 
materials,  the  phase-field  approach  exhibits  a  clear  advantage  in  its  simplicity  in  model  formulation,  ease  of 
numerical  implementation,  and  the  ability  to  explore  essential  interfacial  physics  at  the  interfacial  regions 
etc.  Computing  the  interface  without  explicitly  tracking  the  interface  is  the  most  attractive  numerical  feature 
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of  this  modeling  and  computational  technology.  Since  the  pioneering  work  of  Cahn  and  Hilliard  in  the 
50's  of  the  last  century,  the  Cahn-Hilliard  equation  has  been  the  foundation  for  various  phase  field  models 
[6,  7],  It  arises  naturally  as  a  model  for  multiphase  fluid  mixtures  should  the  entropic  and  mixing  energy 
of  the  mixture  system  be  known.  For  immiscible  binary  fluid  mixtures,  one  commonly  uses  a  labeling  or 
a  phase  variable  (also  known  as  a  volume  fraction  or  an  order  parameter)  <f>  to  distinguish  between  distinct 
fluid  phases.  For  instance  <f)  =  1  indicates  one  fluid  phase  while  4>  =  0  denotes  the  other  fluid  phase  in  an 
immiscible  binary  mixture.  The  interfacial  region  is  tracked  by  0  <  (f)  <  1.  Given  the  historical  reason, 
most  mixing  energies  are  calculated  in  terms  of  the  volume  fraction  instead  of  the  mass  fraction  in  the 
literature  [20,  11].  Consequently,  the  system  free  energy  including  the  entropic  and  mixing  contribution  has 
been  formulated  in  the  volume  fraction  as  well  [20,  11].  We  denote  the  system  free  energy  for  the  material 
system  to  be  modeled  by  77(tj),  V(f),  •  •  • ).  A  transport  equation  for  the  volume  fraction  <f)  along  with  the 
conservation  equation  of  momentum  and  continuity  equation  constitutes  the  essential  paid  of  the  governing 
system  of  equations  for  the  binary  fluid  mixture.  The  volume  fraction  serves  as  an  interval  variable  for  the 
fluid  mixture. 

In  the  literature  on  immiscible  binary  mixtures  of  incompressible  fluids,  one  uses  the  concept  of  chem¬ 
ical  potential  to  formulate  the  transport  equation  for  the  volume  fractions  of  the  fluids  <f)i  and  <[>2 -  In  this 
formulation,  the  material  incompressibility  is  on  the  one  hand  modeled  by  the  continuity  equation 

V  •  v  =  0,  (1.1) 

while  on  the  other  hand,  interpreted  as  the  invariant  property  of  the  sum  of  the  volume  fractions  for  the 
two  fluid  components,  i.e.,  (|)|  +  (>2  =  I  if  we  denote  4>  =  4>t  and  <[>2  =  1  —  ([>-  This  assumption  is  plausible 
and  indeed  consistent  with  the  fluid  compressibility  (1.1)  only  if  the  two  components  are  either  completely 
separated  by  phase  boundaries  when  their  densities  arc  not  equal  or  possibly  mixed  when  the  densities  arc 
identical.  Otherwise,  there  is  a  potential  inconsistency  with  the  usual  conservation  of  mass.  This  inconsis¬ 
tency  has  been  identified  in  [28],  but  ignored  by  many  practitioners  using  phase  field  modeling  technologies. 
We  note  that  this  inconsistency  occurs  only  in  the  mixed  region  of  the  two  incompressible  fluids,  where  the 
incompressible  condition  (1.1)  is  no  longer  valid,  indicating  the  mixture  is  no  longer  incompressible  despite 
that  each  fluid  component  participating  in  mixing  is  incompressible.  This  type  of  fluids  is  referred  to  as 
quasi -compressible  in  [28]. 

This  paper  aims  at  discussing  the  inconsistency  for  binary  mixtures  of  two  incompressible  fluids  of  un¬ 
matched  densities  and  viscosities  and  providing  a  quantitative  assessment  for  the  quasi-compressible  phase 
field  model  that  obeys  the  conservation  of  both  mass  and  volume  against  the  incompressible  one  that  only 
respects  the  volume  conservation.  The  paper  is  organized  as  follows.  First  we  discuss  the  mathematical 
formulation  of  the  phase  field  theory  for  binary  viscous  fluid  mixtures  and  its  various  approximations  and 
their  ramifications.  Secondly,  we  develop  a  new  set  of  numerical  algorithms,  which  enforce  the  mass  con¬ 
servation,  to  solve  the  governing  system  of  fluid  flow  equations.  Thirdly,  we  implement  the  algorithms  using 
spatially  high  order  spectral  methods  and  discuss  the  discrepancies  between  the  ad  hoc  incompressible  phase 
field  model  and  the  quasi-compressible  phase  field  model  in  two  representative  examples. 

2  The  mathematical  model 

We  revisit  the  derivation  of  the  governing  system  of  equations  for  a  binary  mixture  of  incompressible 
viscous  fluids,  which  includes  the  transport  equation  for  a  phase  variable  (the  volume  fraction)  and  the 
conservation  equations  for  mass  and  linear  momentum.  The  conservation  equations  for  a  binary  system  can 
be  formulated  in  two  different  ways:  either  as  a  two  fluid  model  or  a  one  fluid  two  component  model  [1, 
4,  2].  For  nontrivial  fluid  simulations,  the  one  fluid  multi-component  formulation  often  yields  a  convenient 
governing  system  of  equations  and  easy-to-implement  boundary  conditions  for  the  model’s  hydro-dynamical 
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variables.  The  phase  field  theory  falls  naturally  into  the  one  fluid  two  component  formulation  [6].  In  the 
phase  field  formulation,  chemical  reaction  between  the  two  distinct  components  can  take  place  so  that  one 
component  can  be  turned  into  the  other  component.  However,  the  overall  mass  must  be  conserved.  In 
this  paper,  we  will  not  address  the  phase  field  formulation  with  chemical  reactions.  This  topic  deserves  a 
separate  discussion  of  its  own. 

2.1  Governing  equations 

In  a  phase  field  theory,  the  transport  equation  for  the  volume  fraction  of  one  fluid  phase  is  given  by 

<k  +  V  •  (<|)v)  =  V-(XV/u),  (2.1) 

where  v  is  an  average  velocity  to  be  clarified  below,  X  =  X(<|))  is  the  mobility  function,  and  q  =  %  is  the 
chemical  potential  of  the  material  system.  The  mobility  function  X  is  often  taken  as  a  constant  Xo,  but  is 
preferably  a  function  of  (f>  in  the  form: 


A,  =  A.0<1>(1  —  <».  (2.2) 

The  Cahn-Hilliard  equation  with  the  volume  fraction  dependent  mobility  is  called  singular  or  modified 
Cahn-Hilliard  equation.  Often,  it  is  approximated  simply  by  a  constant  value  X  =  Xo  in  studying  phase 
separated,  immiscible  fluids.  The  resultant  equation  is  the  well-known  Cahn-Hilliard  equation. 

The  free  energy  of  the  mixture  system  is  normally  a  function  of  the  labeling  function  of  phase  function 
and  its  higher  order  derivatives  (only  the  first  order  is  included  here  for  brevity): 

F  =  F(<|>,V<|>).  (2.3) 

In  this  paper,  we  consider  the  mixture  of  two  incompressible  fluids  with  constant  mass  density  pi  and  P2, 
respectively.  The  total  density  of  the  mixture  is  then  given  by 

p  =  pi<t>  +  p2(l-<l>).  (2.4) 

We  identify  v  as  the  mass-average  velocity  for  the  mixture.  Then,  the  conservation  equations  for  mass  and 
momentum  arc  given  by 


Pf  +  V  •  (pv)  =  0, 


(2.5) 


(pv),  +  V  •  (pvv)  =  V  •  (x)  -  (|)Vp  +  ¥e, 

where  Fe  is  the  external  force  and  pVp  is  the  ’’elastic  force”  or  the  ’’surface  force”  due  to  the  interfacial 
energy  /(())).  The  surface  force  — (|)Vp  can  be  replaced  by  pV(|)  modulo  a  surface  term  which  is  normally 
zero.  In  light  of  the  transport  equation  for  the  volume  fraction,  we  have 

PV  -  v  =  — (pi  —  p2)(<>r  H- v  -  V<>)  =  —(pi  —  P2)(V  •  (XV /u)  —  <f)V  •  v),  (2.6) 


we  then  derive  from  the  above  and  (2.1)  that 

v  •  V  =  -  [<t>f  +  V-  ( 4>v )]  =  [V  •  (XV/i)] .  (2.7) 

P2  P2 

It  is  apparent  that  the  divergence  free  condition  for  the  mass-average  velocity  field  is  satisfied  only  if  pi  =  p2 
or  V  •  (XV ju)  =  0.  Otherwise,  the  mass  conservation  equation  serves  as  a  constraint  for  the  velocity  field, 
which  determines  the  undetermined  pressure  in  the  constraint  hydrodynamic  theory  for  fluid  mixtures.  We 
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note  that  V  •  (. A,Vp )  is  normally  not  zero  for  a  spatially  inhomogeneous  system.  Hence,  as  long  as  pi  /  p2, 
eq.  (2.7)  serves  as  a  constraint. 

To  close  the  system  of  equations,  we  must  come  up  with  a  constitutive  equation  for  the  stress  tensor  x. 
We  consider  the  mixture  made  up  of  viscous  fluids.  For  viscous  fluids,  the  stress  constitutive  equation  is 

x  =  tc  +  2t|D  +  v/t(D)I,  (2.8) 

where  xc  is  the  constraint  stress  responsible  to  maintain  the  constraint  eq.  (2.7)  without  any  contribution 
to  the  entropy  production,  r|  is  the  shear  viscosity,  v  is  the  volumetric  viscosity,  and  D  is  the  rate  of  strain 
tensor.  The  ratio  between  v  and  r|  depends  on  the  property  of  the  material  and  is  roughly  4.3  for  water  for 
example.  The  viscosity  coefficients  for  the  fluid  mixture  arc  interpolated  through  the  volume  fraction  and 
given  by 


ri  =  nitt>  +  n2(l  —  <t>),v  =  Vi4»  +  v2(l  —  4>),  (2.9) 

where  T|i,2,Vi,2  are  constant  shear  and  volumetric  viscosities  for  fluid  1  and  fluid  2,  respectively. 

To  deal  with  constraint  (2.7),  we  augment  the  free  energy  with  a  term  F  called  the  constraint  response: 

‘f  =  F  +  F.  (2.10) 


Based  on  the  second  law  of  thermodynamics  in  the  form  of  the  Clausius-Duhem  inequality,  the  constraint 
response  does  not  contribute  to  entropy  production,  i.e., 

8F  . 

xc :  D  —  —c|)  =  0,  (2.11) 

where  <j>  =  y  +  V  •  (v(f>)  is  the  material  derivative.  We  rewrite  eq.  (2.7)  as 

I :  D  +  Pl  ~  p2(b  =  0.  (2.12) 

P2 

For  eq.  (2.11)  must  be  valid  for  all  thermodynamic  processes  that  obeys  (2.12),  we  deduce  that 


Tc  =  ~pl  (2.13) 

where  p  is  the  hydrodynamic  pressure,  and  the  free  energy  component  (F)  corresponding  to  the  constraint 
response: 


F  = 


Pl  P2 
P2 


<t >P- 


(2.14) 


If  we  choose  p  =  -g^-,  we  obtain  a  set  of  equations  that  respect  the  conservation  of  mass  and  total  volume: 


<t>r  +  V  •  (<f)v)  =  V  •  (AVp), 

(pv),  +  V  •  (pvv)  =  p[v,  +  v-  Vv]  =  V  •  (2r|D  +  Wr(D)I)  —  Vp  —  ( |)Vp  +  Fe 
=  V  •  (r|Vv)  +  V((r|  +v)V  •  v)  —  Vp  —  (JtVp  +  F^, 

V  •  v  =  [v  •  (AVp)] , 

P2 

=  SfF  =  5F  Pl  —  p2 
P  8(f)  8(f)  '  p2  P ' 


(2.15a) 

(2.15  b) 
(2.15c) 

(2.15d) 
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We  refer  (2.15)  as  the  compressible  model  2  in  this  paper.  On  the  other  hand,  if  we  set  p  =  |^  in  the  system 
of  equations  listed  above,  we  obtain  another  set  of  equations,  which  we  refer  as  the  compressible  model  1. 
With  the  help  of  (2.15c),  the  transport  equation  for  tf»  can  be  recast  into 

§t  +  V  •  (4>v)  =  — — — Vv  (2.16) 

P2-P1 

provided  pi  /  p2- 

The  above  compressible  models  preserve  the  mass  conservation  and  arc  compressible  inside  the  mix¬ 
ing/interfacial  region.  On  the  other  hand,  the  incompressible  model,  in  which  the  mass  average  velocity 
field  is  assumed  solenoidal,  consists  of  the  following  equations: 


<t>r  +  V  -  (4>v)  =  V  •  (A,Vp), 

p[v,  +  v  -  Vv]  =  V  •  (r|Vv))  —  Vp  — <t)Vp  +  Fe, 

V  •  v  =  0, 

8F 

8^' 


(2.17a) 

(2.17b) 

(2.17c) 

(2.17d) 


This  model  assumes  that  the  flow  is  incompressible  everywhere  at  the  expense  of  local  mass  conservation 
inside  the  interfacial/mixing  region. 

For  the  compressible  models,  we  define  the  total  energy  as 

£(t)=  [[^r-  +  f}dx,  (2.18) 

Jq  2 

where  x  is  the  Eulerian  coordinate.  For  compressible  model  1,  the  rate  of  change  in  the  total  energy  is  given 
by 

=  -  [  [^l|Vp||2  +  (T-Tc)  :D]r/x+Pl  P2  f  <]>(v •  Vp+  ^-)dx.  (2.19) 

at  Jq  p2  Jq  at 

To  ensure  positivity  in  the  first  integral  (2r|D  +  Wr(D)I)  :  D  >  0,  r|  >  0,v  +  ^  >  0.  Unfortunately,  we  have 
no  control  over  the  sign  for  the  second  integral.  For  compressible  model  2,  the  rate  of  change  in  the  energy 
for  this  system  of  governing  equations  is  given  by 

? = -//i'v"''2+<^> :  <2-20> 


In  either  models,  the  energy  dissipation  can  not  be  rigorously  established.  On  the  other  hand,  for  the  incom¬ 
pressible  model,  the  energy  is  defined  by 

E{t)=  f[^^+F}dx.  (2.21) 

With  the  help  of  the  divergence-free  condition  in  the  continuity  equation,  Shen  and  Yang  [32,  33]  showed 
an  energy  law  is  available  for  the  incompressible  system: 

d4:  =  -\  [^11  VF||2  +  (X  -  Tc)  :  D }dx  <  0.  (2.22) 

at  Jq, 

One  of  the  traded-offs  in  enforcing  the  mass  conservation  is  losing  the  energy  law.  Therefore,  both  com¬ 
pressible  models  along  with  the  incompressible  mixture  model  in  [33]  arc  valid  approximately.  We  note  that 
an  analogous  phase  field  equation  (to  model  1)  was  also  derived  using  an  energy  argument  in  [28]. 


5 


2.2  Choice  of  free  energy 

The  free  energy  F  can  take  different  form  depending  on  the  applications.  In  this  paper,  we  consider  the 
free  energy  density  in  the  following  form: 

F(^,V^)  =  kBrY[i||V^||2  +  /(^)],  (2.23) 

where  kg  is  the  Boltzmann  constant,  T  is  the  absolute  temperature,  and  y  is  a  parameter  with  the  unit  of  a 
number  density  per  unit  length,  y  is  in  fact  proportional  to  the  product  of  the  number  density  per  unit  volume 
and  the  square  of  the  persistent  length. 

We  first  look  at  the  Ginzburg-Laudau  free  energy  with 

m  =  l-^)2  (2-24) 

for  two  immiscible  fluids,  where  e  >  0  is  a  small  parameter  characterizing  the  hydrophobic  property  between 
the  two  fluids.  Therefore, 


&F  I 

—  =  W[-V2d>  +  ^d>(l  -0>)(1  -2d>)].  (2.25) 

We  also  consider  the  Flory-Huggins  mixing  free  energy  for  two  immiscible  fluids  to  simulate  the  phase 
separation  dynamics.  The  Flory-Huggins  mixing  free  energy  density  is  given  by  (2.23)  with 

fW  =  ~2  Itt""  ln^  ”h  ln(l  (Jt)  +  2C4> ( 1  0)]>  (2-26) 

8  1\\  1\ 2 

where  N\  and  No  are  the  polymerization  indices  for  fluid  1  and  fluid  2  and  %  is  the  mixing  parameter  between 
0  and  2.  If  both  are  viscous  fluids,  we  assume  N\  =  N2  =  1.  In  this  case, 

8F  r  ,  1  ln§  ln(  1  —  d>)  ,  , . , 

—  _  ^B7’y[_v-(|)  +  — (—  -| - — - F)C(1  —  2(f)))]  +  const.  (2.27) 

2.3  Non-dimensionaiization 

We  denote  the  characteristic  time  scale  by  to  and  length  scale  by  Lq.  The  dimensionless  variables  arc 
defined  by 


t  x  vt0  pt?\ 

-,x=—,y=—,p  =  —1. 

to  M)  Lq  p2L() 


(2.28) 


We  will  drop  the ~ on  the  dimensionless  variables  in  the  following.  We  choose  Lq  so  that  the  dimensionless 
length  Ly  =  1.  We  use  Lx  =  1  in  the  following  calculations  simply  for  convenience.  The  dimensionless 
model  parameters  arc  defined  by 


Reu  = 


_  PiLq 

mt 


(i=l,2)A=«jn,  8  =  i 


0  I  1— <l>  1  _  j>  lH 

Re  '  Re2}S  ’  Rev  Re ltV  '  Re2,v  ' 


(2.29) 


Here  Res  and  Rev  denotes  the  Reynolds  number  corresponding  to  the  shear  and  volumetric  stress,  and  A  is 
the  dimensionless  mobility  parameter.  We  set  Ly  =  \  .  C  =  -  =  1  in  this  study  yielding  to  =  jpr,- 
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The  dimensionless  equations  for  the  two  compressible  models  arc  given  by 


+  V  •  (<f>v)  =  V  •  (AVp), 

Bi[\,  +  v  •  Vv]  =  V  •  (j^Vv)  +  V((^  +  ^r)V  •  v)  -  Vp  -  <>Vju, 

<  *  ’  (2.30) 

V-v  =  (l-g)[V-(AVp)], 

k  p  = -AO +  /(0)  (Model  1),  p  =  -A0  +  /(0)  +  (|-l)p(Model2), 

where  /(tj))  arc  given  by  (2.24)  or  (2.26)  with  e  replaced  by  e. 

The  above  system  is  subjected  to  a  set  of  suitable  initial  and  boundary  conditions.  For  example,  if  the 
mixture  is  confined  in  a  domain  £2,  the  boundary  conditions  arc 

9d)  du . 

an  =  5- 0Q  =  0,  vaQ  =  0,  (2.31) 

an  an 

where  n  is  the  outward  normal. 


3  Numerical  schemes 

Shen  and  Yang  have  developed  and  studied  efficient  numerical  schemes  for  the  incompressible  phase 
field  model  in  great  detail  in  [32,  33].  In  this  section,  we  focus  on  extending  their  results  to  the  compressible 
models  given  for  example  by  (2.15)  for  the  binary  fluid  mixture.  Since  the  two  compressible  models  do 
not  admit  a  dissipative  energy  law,  it  is  not  possible  to  construct  energy  stable  schemes  for  them.  Instead, 
we  shall  construct  efficient  and  easy  to  implement  numerical  schemes  to  accommodates  the  non-vanishing 
divergence  velocity  field  so  as  to  preserve  the  mass  conservation. 


3.1  Discretization  in  time 

To  simplify  the  presentation,  we  shall  present  only  first-order  schemes.  In  what  follows,  the  superscript 
n  denotes  the  time  level  and  At  is  the  time  step  size. 


Scheme  based  on  a  modified  projection: 


1.  Solve  (C+1,p"+1)  from: 


AM-Tl  _  /K» 

-  ^  +V-  (<f  v”)  =  V  •  (AVp',+1), 

p”+1  =  — AC+1  +  /(C)  +  p(C+1  -  C) 

where  S  is  a  computational  parameter  and  e  is  the  parameter  in  the  free  energy.  The  last  term  is  added 
to  stabilize  the  scheme  to  allow  larger  step  sizes.  Its  role  is  to  damp  the  high  frequency  or  short  waves 
in  the  numerical  simulation. 


0(|)''+1 

dn 

dn 


las  —  0, 


Ian 


=  0, 


(3.1) 


2.  Denote 

fli"+1=C+1  —  +  (1-C+1),  Re"+1  =<bn+i Ret  s  + (l- <bn+1)Re2s, 

P2  ’  ’  (3.2) 

Ren+l  =  C+1fo1)V  +  (l  -C+1)^2,v; 
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Solve  v'1+1  from: 


5/"+1  ( - — - f  v"-Vv"  j  —  V  •  (/?e"+1Vv"+1)  —  V((/?e"+1  +  /te"+1)V  •  v")  +  Vpn  =  -§n+lVpn+1, 

v'1+1|9r2  =  0; 

(3.3) 


3.  Set  co  =  1  and  solve  pn+l  —  p"  from 

1 


-V- 


n+ 1  ^  j  ^  P2  Pi 


^  V(p',+1  -pn)  =  —  {  c0 
fltn+l  U  F  '  At  l  p2 


[V  •  (AV//+1 )]  —  V  •  v,,+1 1 , 


d(pn+l  —  pn) 
d  n 


(3.4) 


0£2  —  0, 


where  £2  is  the  domain  occupies  by  the  fluid  mixture. 
4.  Finally,  update 


v"+i  =  v"+1  - 


At 

Bi"  1 


V(p'!+1-p”), 


(3.5) 


and  then  goto  the  next  step. 

Remarks: 

•  S  =  0(1)  is  a  stabilizing  computational  parameter.  We  use  S  =  2  in  all  the  simulations  presented  in 
this  paper. 

•  Setting  co  =  0  in  (3.4),  we  get  the  scheme  for  divergence-free  velocity  field  (V  •  v  =  0). 

•  ri  =  <(>ri i  +  (1  —  <|))ri2  and  v  =  cf)V i  +  (1  —  (f))v2  are  the  interpolated  effective  viscosity  coefficients. 

•  A  second-order  scheme  can  be  constructed  as  well. 


Notice  that  (3.4)-(3.5)  represents  a  modified  pressure-correction  projection  method.  One  can  easily 
verify  from  (3.4)-(3.5)  that  v',+1  and  pn+1  satisfy 

V  •  v',+1  =  (1-— )[V-(A  Vp'1+1)],  (3.6) 

P2 


which  ensures  the  mass  conservation.  However,  the  step  (3.4)  in  the  above  involves  solving  an  elliptic  equa¬ 
tion  with  as  the  variable  coefficient.  When  is  large,  this  step  may  become  very  costly.  So  we  propose 
the  following  scheme  based  on  the  pressure-stabilization  technique  which  only  requires  solving  a  pressure 
Poisson  equation.  The  price  we  pay  for  this  simplicity  is  that  (3.6)  will  only  be  satisfied  approximately. 
This  strategy  has  been  proven  effective  in  the  numerical  solution  of  the  incompressible  field  phase  model 
[32,  33], 

Scheme  based  on  a  pressure-stabilization  method: 

1.  Solve  (<|)''+1,//+1)  from: 


<t)"+1-<t)” 


+  V.(<fv")  =  V-(AVq',+  1), 


At 


S 


3(|) 


n+ 1 


dn 


-|an  =  0, 


dp. 


=  -Ar+1 +/(o + ^ (c+1  -<n,  = o. 


(3.7) 


2.  Denote 


Bin+l  =C+1  — +  (1— f'+1),  Re"+l  =f+lRels  +  {l-f+1)Re2s, 

P2  ’  ’  (3.8) 

Renv+i =r+1^1)V+(i-r+1)^2,v; 


Solve  v',+1  from: 

Bin  1  1  ( - — - +  v"  •  Vv"  J  -  V  •  (Ren+lW+1)  -  V((Rens+1  +  tfe”+1)V  •  v”)  +  Vpn  =  -(|)',+1  Vp'1+1 , 

v'1+1|an  =  0; 

(3.9) 


3.  Set  Co  =  1  and  solve  pn+]  —  p”: 


d(pn+l-p ») 
dn 


At 


las  —  05 


-  A(p"+1  ^  i  c0P2  Pl[V  •  (AV/+1)]  -  V  •  v',+1 


P2 


(3.10) 


where  pm/>,  =  min(pi,p2).  Go  to  the  next  step. 


We  observe  that  the  step  (3.7)  is  a  system  of  two  second-order  equations  with  constant  coefficients,  the 
step  (3.9)  is  an  elliptic  equation  with  variable  coefficients  and  the  step  (3.10)  is  just  a  Poisson  equation. 
Hence,  the  above  scheme  is  easy  to  implement  and  very  efficient. 

As  in  an  usual  pressure-stabilization  method  [33]  where  the  divergence-free  condition  is  satisfied  ap¬ 
proximately,  it  is  clear  from  (3.10)  that  v',+1  and  p" + 1  from  the  above  scheme  only  satisfy  the  internal 
constraint  (2.7)  approximately  with  a  residue  of  order  0(At2).  Therefore,  the  mass  is  conserved  up  to  a 
controllable  error  of  order  0(At2),  indepedent  of  the  interfacial  width  e. 


3.2  Discretization  in  space 

The  spatial  discretization  can  be  done  in  either  a  spectral  method  or  a  finite  element  method  or  a  finite 
difference  method.  However,  the  spatial  resolution  needs  to  be  fine  enough  to  resolve  the  interfacial  layer. 
We  shall  use  the  high  resolution  spectral  method  which  requires  a  significantly  less  number  of  unknowns 
inside  the  interface  as  compared  with  a  lower-order  method. 

We  focus  in  this  paper  on  two-dimensional  fluid  flows  in  both  drop  dynamics  as  well  as  mixing  dynamics 
of  immiscible  binary  fluids,  and  define  the  computational  domain  as  £2  =  [0,LV]  x  [0,LV]  with  the  periodic 
boundary  condition  in  the  x-direction.  In  the  v-direction,  the  boundary  conditions  arc: 

\(x,y,t)  =y(x  +  Lx,y,t),<$>(x,y,t)  =  §(x  +  Lx,y,t), 

(3.11) 

v(x,0  ,t)  =  V0  ,\(X,Ly,t)  =  Vl  ,  (J)-y  (x,  0 ,  f )  =  tyy(x,Ly,t)  =  ^>yyy  (x,  0  ,t )  =  tyyyy(x,Ly,t)  =  0. 

The  boundary  conditions  of  (|)  at  y  =  0,Ly  arc  interpreted  as  the  flux  boundary  conditions. 

We  shall  use  the  Fourier  expansion  in  the  x-direction  and  the  Legendre-Galerkin  method  in  the  y- 
direction. 
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Table  1 :  Parameter  values 


Parameter 

N 

M 

8 

Rei,s 

Rei,s 

e 

Rei,v 

R&2,v 

(3i 

P2 

A 

Value 

256 

256 

0 

1 

1  x  10~2 

0.02 

4.3  x  ReitS 

4.3  x  Re 2,v 

1 

50 

1  x  10~b 

4  Numerical  results  and  discussions 

We  investigate  predictive  drop  dynamics  and  phase  separation  dynamics  computed  using  the  two  distinct 
classes  of  models  with  a  focus  on  the  comparison  between  the  phases  and  the  phase  boundaries  of  the 
mixture.  We  tabulate  the  dimensionless  parameters  used  in  the  simulations  in  Table  1,  where  N  and  M 
denote  the  number  of  grid  points  in  x  and  y  direction,  respectively.  These  arc  chosen  based  on  our  previous 
experience  with  the  two-phase  fluid  [32,  33]. 

We  first  consider  the  drop  dynamics  of  fluid  1  immersed  in  fluid  2  and  denote  the  volume  fraction  of  fluid 

1  as  <f).  For  presentation  purposes,  we  relabel  the  models  as  follows  in  the  figures:  Model  1:  incompressible 
model.  Model  2:  compressible  model  1,  Model  3:  compressible  model  2. 

4.1  Drop  dynamics 

We  first  simulate  a  lighter  fluid  (fluid  1 .  <|)  =  1 )  drop  immersed  in  a  heavier  fluid  (fluid  2,  <|)  =  0).  The 
density  ratio  we  choose  for  this  numerical  example  is  pi  :  P2  =  1  :  50  and  viscosity  ratio  1  :  100.  In  this 
setting,  the  lighter  drop  will  rise  in  the  fluid  channel  (computational  domain).  The  simulated  results  using 
the  compressible  models  and  the  incompressible  model  for  mixtures  agree  with  each  other  qualitatively. 
The  velocity  components,  pressure  and  the  drop  profiles  obtained  using  the  three  phase  field  models  arc 
shown  in  Figures  1-3,  respectively.  In  the  simulations,  the  lighter  fluid  drop  rises;  the  rising  drop  pushes 
the  fluid  in  the  front  aside  and  pushes  the  fluid  downward  on  the  side  of  the  fluid  domain.  The  horizontal 
and  the  vertical  velocity  component  are  plotted  in  Figures  1  and  2,  respectively.  The  pressure  around  the 
drop  remains  low,  which  is  shown  in  Figure  3.  The  drop  shapes  obtained  using  the  three  distinct  models  arc 
contrasted  at  a  selected  time  t  =  3.6  in  Figure  4  along  with  the  deviations  between  the  velocity  components 
of  each  pair  of  models.  The  predictions  from  model  1  and  2  are  close  relative  to  that  from  model  3.  The 
deviations  in  general  fall  into  the  range  of  (9(10  2).  The  velocity  field  superimposed  by  the  drop  profile  is 
shown  in  Figure  5,  where  a  pair  of  vortices  arc  shown  explicitly. 

We  then  repeat  the  simulation  with  a  heavier  fluid  drop  sediments  in  a  lighter  fluid.  The  density  ratio  is 
reversed  to  pi  :  P2  =  50  :  1  and  the  viscosity  ratio  is  reversed  to  100  :  1.  The  behavior  described  above  for 
the  rising  drop  reverses.  This  time,  the  predictions  between  model  1  and  model  2  and  those  between  model 

2  and  model  3  are  qualitatively  the  same;  model  3  predicts  the  fastest  drop  sedimentation  among  all  three. 
To  save  space,  we  suppress  the  demonstration  of  the  numerical  results  pertinent  to  this  simulation. 

In  summary,  the  model  predictions  in  both  drop  rising  and  drop  sedimentation  agree  qualitatively.  In 
these  cases,  the  interfacial  layer  between  the  two  immiscible  fluids  arc  thin  and  the  volume  fraction  of 
the  fluid  involved  in  the  mixing/interfacial  zone  is  small.  Consequently,  the  deviation  among  the  model 
predictions  arc  small.  We  anticipate  this  scenario  will  change  as  the  mixing/interfacial  layer  gets  larger  and 
the  volume  fraction  of  the  fluid  involved  in  mixing  becomes  significant.  We  next  examine  an  example  of 
fluid  mixing/demixing  where  the  mixing  zone  is  significantly  larger. 

4.2  Phase  separation  dynamics  of  immiscible  binary  fluids 

Figures  6-9  depict  the  phase  portrait  of  the  mixture  during  phase  separation  and  the  corresponding 
velocity  components  as  well  as  the  pressure  field  at  selected  time.  In  Figure  6,  the  value  of  the  volume 
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fraction  cf)  is  plotted  as  a  color  map.  The  compressible  models  (model  2  and  model  3)  give  well  separated 
islands  while  the  incompressible  model  (model  1)  predicts  only  slightly  modified  phase  landscape.  Figures  7 
and  8  supports  this  with  a  much  elevated  velocity  field  in  the  compressible  models  than  in  the  incompressible 
model.  Moreover,  the  flow  pattern  is  drastically  different  between  the  predictions  obtained  from  different 
classes  of  models.  Figure  9  portraits  the  pressure  field,  which  correlates  well  with  the  phase  portrait  of 
the  mixture  given  by  the  level  sets  of  the  volume  fraction  (f>.  The  difference  between  the  two  classes  of 
models  are  significant  in  this  numerical  example.  The  drastic  difference  between  the  model  predictions  is 
an  amplification  of  the  difference  in  the  fundamental  physical  mechanism  on  mass  conservation  in  a  much 
larger  mixing  zone  in  contrast  to  the  previous  drop  dynamics. 

If  these  examples  show  the  behavior  of  the  transient  solution,  the  next  set  of  figures  (Figures  10-13) 
portrait  the  solutions  up  to  nearly  quasi-static  states.  The  phase  behavior  predicted  by  the  incompressible 
model  (Model  1)  is  distinct  quantitatively  from  those  by  the  compressible  models  (Figure  10).  The  prediction 
on  the  velocity  field  and  the  pressure  field  made  by  the  incompressible  model  and  by  the  compressible  ones 
are  completely  different.  Whereas,  the  difference  between  the  compressible  model  predictions  is  minimal. 

If  we  were  to  impose  the  constraint  on  the  conservation  of  the  total  volume  of  each  separate  phase, 
the  predictions  from  the  compressible  models  should  be  more  credible  since  they  also  conserves  the  mass, 
which  is  fundamentally  important. 

5  Conclusion 

A  pair  of  phase  field  models  that  conserve  mass,  momentum  and  total  volume  for  each  individual  phase 
of  immiscible  binary  fluid  mixtures  are  formulated.  The  mass-average  velocity  becomes  non-solenoidal 
when  the  density  ratio  between  the  two  fluids  is  not  unity.  Consequentially,  the  new  phase  field  theories 
are  compressible  although  a  global  volume  conservation  for  each  phase  can  be  maintained  over  the  entire 
material  volume.  The  commonly  used  phase  field  model  for  binary  fluid  mixtures  is  the  model  we  refer  to 
as  the  incompressible  model  in  this  paper,  in  which  the  continuity  equation  is  approximated  by  a  divergence 
free  condition;  the  resulting  theory  preserves  the  material  volume  but  not  the  mass. 

The  deviation  between  the  predictions  by  the  compressible  models  and  the  incompressible  one  depend 
on  the  size  of  the  mixing  zone.  When  the  size  of  the  mixing  zone  is  small  compared  to  the  entire  fluid 
domain,  the  model  predictions  agree  qualitatively.  However,  when  the  mixing  zone  is  large,  the  two  classes 
of  models  describe  two  quite  different  dynamics  (in  both  transient  and  quasi-steady  state).  One  numerical 
example  on  a  drop  dynamics  of  one  fluid  drop  immersed  in  another  immiscible  fluid  matrix  and  the  other 
on  the  phase  separation  of  immiscible  binary  fluid  mixtures  arc  carried  out  to  illustrate  the  point.  From  the 
hydrodynamics  point  of  view,  it  is  apparent  that  the  fundamental  conservation  laws  of  fluids  must  be  obeyed. 
Therefore,  the  mass  conservation  should  be  respected  in  any  faithful  simulations  employing  the  phase  field 
formulation  when  the  mixing  zone  is  large.  The  predictions  made  by  the  compressible  models  arc  consistent 
in  the  two  numerical  examples,  and  therefore  arc  credible  regardless  of  the  size  of  mixing  zone. 
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Figure  1:  The  horizontal  velocity  component  in  the  case  of  a  lighter  drop  immersed  in  a  heavier  fluid  with 
density  ratio  pi  :  p2  =  1  :  50  at  time  t  =  0.6, 1.2, 1.8, 2.4, 3, 3.6.  (a)  Ml  (Model  1):  (b).  M2  (Model  2):  (c). 
M3  (Model  3).  Each  figure  is  superimposed  by  the  shape  of  file  drop,  i.e.  the  zero  contour  curve  of  the 
phase-field  function  <|),  at  the  corresponding  time. 
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(a) 


Figure  2:  The  vertical  velocity  component  in  the  case  of  a  lighter  drop  immersed  in  a  heavier  fluid  with 
density  ratio  pi  :  P2  =  1  :  50  at  t  =  0.6, 1.2, 1.8, 2.4, 3, 3.6.  (a)  Ml,  (b).  M2,  (c).  M3.  Each  figure  is 
superimposed  by  the  shape  of  the  drop,  i.e.  the  zero  contour  curve  of  the  phase-field  function  <|>,  at  the 
corresponding  time. 
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Figure  3:  The  pressure  field  in  the  case  of  a  fighter  drop  immersed  in  a  heavier  fluid  with  density  ratio 
pi  :  p2  =  1  :  50  at  t  =  0.6, 1.2, 1.8, 2.4, 3, 3.6.  (a)  Ml,  (b).  M2,  (c).  M3.  Each  figure  is  superimposed  by  the 
shape  of  the  chop,  i.e.  the  zero  contour  curve  of  the  phase-field  function  <(),  at  the  corresponding  time. 
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Figure  4:  The  velocity  difference  between  pairs  of  die  phase  field  models  investigated,  (a)  h(M1)  —  u(M2), 
(b)  u(M\)  -  u (M3),  (c)  u(M2)  -  u(M3),  (d)  v(Ml)  -  v(M2),  (e)  v(Ml)  -  v(M3),  (f)  v(M2)  -  v(M3),  (g) 
Comparisons  of  the  zero  contour  curves  of  (j>  for  all  three  models.  The  results  are  based  on  the  solutions  at 
t  =  3.6. 
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Figure  5:  Tire  velocity  field  superimposed  by  the  shape  of  the  drop  at  t  =  3.6  for  die  three  models,  (a)  Ml 
(b)  Ml ,  (c)  M3. 
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Figure  6:  Tire  phase  portrait  of  the  binary  fluid  mixing  with  the  Flory-Huggins  energy  for  all  three  models 
at  t  =  0.1, 0.2, 0.3, 0.4, 0.5.  (a)  Ml,  (b)  M2,  (c)  M3. 
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Figure  7:  The  velocity  component  u  at  t  =  0.1, 0.2, 0.3, 0.4, 0.5  in  the  case  of  binary  fluid  mixing  for  all  three 
models,  (a)  Ml,  (b)  M2,  (c)  M3. 
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Figure  8:  The  velocity  component  v  at  t  =  0.1, 0.2, 0.3, 0.4, 0.5  in  the  case  of  binary  fluid  mixing  for  all  three 
models,  (a)  Ml,  (b)  M2,  (c)  M3. 
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Figure  9:  The  hydrodynamic  pressure  p  at  t  =  0.1, 0.2, 0.3, 0.4, 0.5  in  the  case  of  binary  fluid  mixing  for  all 
three  models,  (a)  Ml,  (b)  M2,  (c)  M3. 
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Figure  10:  The  phase  portrait  of  the  binary  fluid  mixing  with  the  Floiy-Huggins  energy  for  all  three  models 
at  t  =  0.1,0.2,03,0.4, 15.  (a)  Ml,  (b)  M2,  (c)  M3. 
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Figure  11:  The  velocity  component  u  at  t  =  0.1, 0.2, 0.3, 0.4, 15  in  the  case  of  binary  fluid  mixing  for  all 
three  models,  (a)  Ml,  (b)  M2,  (c)  M3.  The  solutions  nearly  reach  quasi-steady  state  at  t  =  15. 
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Figure  12:  The  velocity  component  v  at  t  =  0.1, 0.2, 0.3, 0.4, 15  in  the  case  of  binary  fluid  mixing  for  all 
three  models,  (a)  Ml,  (b)  M2,  (c)  M3.  The  solutions  nearly  reach  quasi-steady  state  at  t  =  15. 
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Figure  13:  The  hydrodynamic  pressure  p  att  =  0.1, 0.2, 0.3, 0.4, 15  in  the  case  of  binaiy  fluid  mixing  for  all 
three  models,  (a)  Ml,  (b)  M2,  (c)  M3.  The  solutions  nearly  reach  quasi-steady  state  at  t  =  15. 
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